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We study the complete cosmological evolution of dark matter (DM) density fluctuations for DM 
particles that decoupled being ultrarelativistic during the radiation dominated era which is the 
case of keV scale warm DM (WDM). The new framework presented here can be applied to other 
types of DM and in particular we extend it to cold DM (CDM). The coUisionless and linearized 
Boltzmann-Vlasov equations (B-V) for WDM and neutrinos in the presence of photons and cou- 
^ ' pled to the linearized Einstein equations are studied in detail in the presence of anisotropic stress 

with the Newtonian potential generically different from the spatial curvature perturbations. We 
recast this full system of B-V equations for DM and neutrinos into a system of coupled Volterra 
integral equations. These Volterra-type equations are valid both in the radiation dominated (RD) 
and matter dominated (MD) eras during which the WDM particles are ultrarelativistic and then 
nonrelativistic. This generalizes the so-called Gilbert integral equation only valid for nonrelativistic 
particles in the MD era. We succeed to reduce the system of four Volterra integral equations for the 
density and anisotropic stress fluctuations of DM and neutrinos into a system of only two coupled 
■ Volterra equations. The kernels and inhomogeneities in these equations are explicitly given func- 

tions. Combining the Boltzmann-Vlasov equations and the linearized Einstein equations constrain 
the initial conditions on the distribution functions and gravitational potentials. In the absence of 
neutrinos the anisotropic stress vanishes and the Volterra-type equations reduce to a single integral 
^ ■ equation. These Volterra integral equations provide a useful and precise framework to compute the 

i—^ ' primordial WDM fluctuations over a wide range of scales including small scales up to ~ 1/5 kpc. 
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I. INTRODUCTION AND SUMMARY OF RESULTS 

The evolution of the dark matter (DM) density fluctuations since the DM decoupling till today is a basic problem in 
cosmology. This problem has been extensively treated in the literature for particles decoupling being nonrelativistic 
(Cold DM, CDM) ;i-5]. 

Particles decoupling ultrarelativistically in the radiation dominated era (Warm DM, WDM) were proposed as DM 
candidates year s ago ■ Such WDM particles with mass in the keV scale become the subject of a renewed interest 
in recent years |10l - ll6l |. 

In this paper, we study the evolution of DM density fluctuations for particles that decoupled being ultrarelativistic 
during the radiation dominated era. (Ref. [i3 has recently considered this issue). 

The expansion of the Universe dilutes matter in the early universe and particle decoupling happens when the particle 
collisions become sufficiently rare and can be neglected. Therefore, and it is well known, the particle distribution 
generically freezes out at decoupling. This happens irrespective of whether the particles are in or out of thermal 
equilibrium (see ref. section 2 of ref.fl^ and ref. [2l|) 

The treatment of the cosmology density fluctuations presented here and in the companion paper ref. ^35^ is valid 
for generic freezed out distribution functions, whether at thermal equilibrium or out of thermal equilibrium and holds 
irrespective of the particular DM particle model. 

The linearized Boltzmann-Vlasov (B-V) equation provides an appropriate framework to follow the evolution of the 
primordial fluctuations since the DM decoupling till today. The linearized B-V equation turns to be particularly 
difficult to solve since it is in general a partial differential equation on a distribution function which depends on seven 
variables. Two strategies have been used to solve the linearized B-V equation. One method consists in expanding the 
distribution function on Legendre polynomials transforming the B-V equation into an infinite hierarchy of coupled 
ordinary differential equations (ODE) [H-ll]. Another approach to the linearized B-V equation integrates the distri- 
bution function over the particle momenta and recast the linearized B-V equation into a linear integral equation of 
Volterra type [20l - |22l . [26j . In the case of nonrelativistic particles in a matter dominated universe this leads 

to the so called Gilbert equation [l^. This approach leads to linear integral equations of Volterra type while the 
Legendre polynomials expansion produces an infinite hierarchy of coupled ODE. The Volterra type integral equation 
exhibits a long-range memory of the gravitational interaction ,26,] . However, the memory of the RD era turns out to 
fade out substantially in the MD era. 

In this paper we derive a system of integral equations of Volterra type valid for relativistic as well as for non- 
relativistic particles propagating in the radiation and matter dominated eras. For warm dark matter and neutrinos 
we obtain a pair of coupled Volterra integral equations for the density fluctuations and the anisotropic stress. 

We start by writing down the collisionless Boltzmann-Vlasov equation in a spatially flat FRW spacetime with 
adiabatic fluctuations in the conformal gauge. The distribution function /rf,„(7/, (f, x) of the DM particles after their 
decoupling and to linear order in the fluctuations can be written as 

UniV, q, = Nam gam /o*"(<Z) + ft"\V, q, X) = Ndm /o*"(g) 9drn [l + ^ dmiV, 9, ^)] , (1-1) 

where /□'"(q) is the homogeneous and isotropic zeroth order distribution at decoupling, gdm is the number of internal 
degrees of freedom of the DM particle and Ndm is a normalization factor, t] is the conformal time, q and x stand 
for the particle momentum and position, respectively. We use the superscript tilde in conflguration space as ^(x) to 
indicate the Fourier transform of the momentum space function ^(fc). The superscript hat stands for dimensionful 
functions as /o™(g) whose dimensionless counterpart /o™(Q) does not bear a hat. 

Neutrinos are analogously described by a distribution function /i/(77, q, x) 

Uif], q, x) = N,{r,) f^iq) + f^{r^, q, f) = N,{r,) /^(q) [l -|- §,(r,, q, x)] , (1.2) 

where /g (?) stands for the zeroth order Fermi-Dirac distribution function for neutrinos after decoupling, g^ is the 
number of neutrino internal degrees of freedom and N^{r]) is a normalization factor. 
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We obtain as the collisionless B-V equation for DM including linear terms in the fluctuations 



—p. h — Oi^dm H 77\ — 

orj h omq 



dm 




dct^ E ^ 

-T. n QtOlV 

orj 



= 



(1.3) 



The neutrino distribution function obeys the massless version of eq. (|1.3l) 

Hi di tp 



+ Ui di + 



dr] 



din g 



= 



(1.4) 



where tp is the Newtonian potential and corresponds to the spatial curvature perturbation. 

The B-V equations (ll.3p - (|1.4p are coupled to the linearized Einstein equations for the gravitational potentials ip 
and (p. After Fourier transforming, the linearized Einstein equations read 



3 h{i^) + fc2 0(r/, ^) + 3 h^i-n) ^j{r], k) = -4tt G 



Adm{v,k)+A^{r),k) 



+ 4a2(y7) p^{r]) 60(77, fc) 



(7(77, fc) = 0(77, fc) — 7/^(77, fc) 



47r G 
fc^ 0-^(77) 



^dmiV, k) + S^(77, fc) = <7dm{rii k) + (Ti,(?7, fc) , 



, (1.5) 
(1.6) 



where Pj^r]) is the photon energy density, 60(77, fc) the photon temperature fluctuations integrated over the q directions, 
ct(77, fc) is the anisotropic stress perturbation and 



Ad„(77, fc) ^ / ^ i?(77, q) /*"(77, g, fc) , J^dmiv, k)^-2 [ ^ ^^^^ P2 (fc • q) 7^(77, g, fc) 



(2 7r)3 



(2^)3 i?(77,g) 



1 da 



hii]) — - ^ , E{ri,q) = J ra^ o?{ri) + cp- 
a drj 



(1.7) 



P2{x) the Legendre polynomial of order two. Equations analogous to eq. (|1.7p hold for neutrinos with the index d„ 
replaced by i, and E{r/,q) replaced by q. 

The customary DM density contrast 6(7], fc) is connected with AdmiVi k) by 



S{r/,k) 



where pdm is the average DM density today. 



Adrn{ri, fc) 



Pdm {P'cq 



3200 



(1.8) 



We start this paper by deriving the collisionless Boltzmann-Vlasov equation for DM particles which decoupled 
being ultrarelativistic (UR) and become non-relativistic in the radiation dominated era. This treatment is general 
and applies to any DM particle candidate decoupling being UR during the RD era. In particular, it is appropriated 
for keV scale WDM particles which become non-relativistic by redshift z ^ 5 x 10^. Furthermore, we generalize the 
whole treatment to particles that decouple being non-relativistic as CDM. 



Combining the linear and collisionless Boltzmann-Vlasov equations (|1.3l) - (|1.4p with the linearized Einstein equations 
p.5p - (|1.6p at initial times strongly constrain the initial conditions on the distribution functions and the gravitational 
potentials. The initial conditions are efficiently investigated expanding the distribution functions in powers of 77 and 
i q ■ k rj. Our analysis includes the initial conditions for DM, neutrinos and photons. This analysis is valid both 
for DM that decouples being ultrarelativistic and nonrelativistic (as CDM). We show in this framework that the 
fc dependence factorizes out in the initial distribution functions ^'^^(O, (f, fc) and ^jy(0, g, fc) as well as in the initial 
densities Acj„i(0, fc), A^(0, fc) and anisotropic stresses ffdmiO, fc), a^{0, fc). The dependence on the directions of fc stays 
factorized for all times considerably simplifying the dynamical evolution. 

The primordial inflationary fluctuations [H, [l^l determine the initial gravitational potential ^j{0, fc). 7/;(0, fc) is given 
by the product of a fc dependent amplitude proportional to times a Gaussian random field with unit variance 

that depends on the fc-direction, 77^ being the scalar primordial index. 



We derive from the linearized Boltzmann-Vlasov equation (|1.3|) a system of four linear integral equations of Volterra 
type for the density fluctuations Adm{il, fc), ^,^(77, fc) and the anisotropic stress fluctuations adm{ih k), <^u{fl, fc) valid 
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both for ultrarelativistic and non-relativistic particles in the RD and MD eras. This is a generalization of Gilbert's 
equation. Gilbert's equation is only valid for non-relativistic particles in a matter dominated universe [T9| . The 
remarkable fact in these new Volterra integral equations is that the density and anisotropic stress fluctuations obey 
a closed system of integral equations. Although the B-V equation is an equation on functions of rj, k and q with 
coefficients depending on 77, k and q, integrating the distribution functions on q with appropriated weights, the 
density and anisotropic stress fluctuations obey a closed system of integral equations. Namely, no extra information 
on the q dependence of the distribution functions is needed, which is a truly remarkable fact. 

In summary, the pair of partial differential Boltzmann-Vlasov equations in seven variables eqs. (|1.3p and (|1.4p 
becomes a system of four Volterra linear integral equations on Adm(77, fc), Edm(?7, fc), ^1,(77, fc) and Y,^{ri,k). In 
addition, because we deal with linear fluctuations evolving on an homogeneous and isotropic cosmology, the Volterra 
kernel turns to be isotropic, independent of the k directions. As stated above, the k dependence factorizes out and 
we arrive to a final system of two Volterra integral equations in two variables: the modulus k and the time that we 
choose to be the scale factor. 

We have thus considerably simplified the original problem: we reduce a pair of partial differential B-V equations 
on seven variables 77, q, x into a pair of Volterra integral equations on two variables: k. 

It is convenient to define dimensionless variables as 




lfs = — 



2 Td / /f " 



where I jg stands for the free-streaming length [l^, [l^ , Td is the comoving DM decoupling temperature and if"^ 
is the dimensionless square velocity dispersion given by 

= / Q" /o™(Q) dQ , while /o''™(Q) is normalized by I^"' = 1 . (1.9) 
Jo 

Q is the dimensionless momentum Q = q/Td whose typical values are of order one. We choose as time variable 

y EE a{ri)/aeq ^ 3200 a{rj) . 

A relevant dimensionless rate emerges: the ratio between the DM particle mass m and the decoupling temperature 
at equilibration, 



_ maeq _ m r gd \ 

= ^ - 4900 — ^— j 



gd being the effective number of UR degrees of freedom at the DM decoupling. Therefore, S^dm is a large number 
provided the DM is non-relativistic at equilibration. For m in the keV scale we have S^dm ^ 5000. 

DM particles and the lightest neutrino become non-relativistic by a redshift 

ztrans + 1 = ^ ^ 1.57 X lo'^ ( T^l " ^r DM particles , Zf^^^^ = 34 ^ "^"^^ for the lightest neutrino . 

Id keV VlOO/ 0.05 eV 

(1-10) 

Ztrans denoting the transition redshift from ultrarelativistic regime to the nonrelativistic regime of the DM particles. 

Since WDM decouples being ultrarelativistic it contributes to radiation for large redshifts z > Ztj-ans- However, 
WDM turns to produce a small relative correction of the order 1 / £^dm to the photons -I- neutrino density. We find a 
little slow down of the order l/^dm in the expansion of the universe when the WDM becomes non-relativistic around 

We obtain a pair of coupled Volterra equations for the functions A(y, a) and a{y, a) defined as follows 



A,„(y,a) + ^A.(y,a) 

Srfm -'3 



. o-(y,a) = 0(?/,a)-l , (t){y,a) = ' ^ , A(0,a) = l, 

7p{0,k) 

(1.11) 

where we choose to factor out the initial gravitational potential V'(0, fc), Rv{y) is the neutrino fraction of the average 
energy density, 

7^™ - TO - - IK - 

h = -^+RM , Ad,„(2/,a) = Admiv,k) , A,(y,a) = A,(r,, fc) , 

?rfm pdrn Td 7p(0, k) Pr V'(0, fc) Ryijj) 



5 



^dmiVik) being given by eq. (|1.7|) . Notice that the DM contribution to A{y,a) is suppressed by a factor 1/^^ 
1/5000. The growth of the DM fluctuations in the MD era largely overcomes this suppression. 



Expanding the Boltzmann-Vlasov equations (jl.3p and (II. 4p in powers of rj and i q ■ k rj as remarked above, we 
obtain the initial gravitational potentials related by 



mk) = 



1 



^(o,fc) 



1 



i?.(0) 



mk) 



Notice above the small 0{l/£,dm) correction in I^. 

The final pair of dimensionless Volterra integral equations take the form 



A{y, a) = C(2/, a) + B^iy) 4>{y, a) + / dy' [G„(j/, y') 4>iy' , a) + G^(y, y') a{y' , a)] , (1.12) 

Jo 

a(y,a) = C"^(y,a)+ r dy' [l^iy,y') aiy'^a) + Uy,y') ^{y' ,a)] , (1.13) 
Jo 

with initial conditions A(0, a) = 1 , (7(0, a) = | . This pair of Volterra equations is coupled with the linearized 
Einstein equations (ll.5l) - (|1.6p . 

The kernels and the inhomogeneous terms in eqs. (ll.l2l) - (|1.13l) are given explicitly by eqs. (|4.7l) - (|4.10l) . (j4.15p - (|4.17p 
and (|4.29p - (|4.33p . The arguments of these functions contain the dimensionless free-streaming distance l{y,Q), 



l{y.Q) 



dy' 



{Q/i 



din ) 



(1.14) 



The coupled Volterra integral equations (|1.12p - (|1.13p are easily amenable to a numerical treatment. 

When the anisotropic stress (j{y,a) is negligible, eas. (ll.l2l) - (|1.13l) reduce to a single Volterra integral equation for 
the DM density fluctuations Adm{y,C() . We find the solution of this single Volterra equation for a broad range of 
wavenumbers 0.1 /Mpc < fc < 1/5 kpc and analyze the transfert function and density contrast in the accompanying 
paper ref. [3^ . 

The framework derived here reducing the full evolution of the primordial cosmological fluctuations to a pair Volterra 
integral equations is general for any type of DM and provides in particular, in the nonrelativistic limit in the MD era 
the so-called Gilbert equations. 

In summary, the Volterra integral equations obtained here provide a useful and precise method to compute the 
primordial DM fluctuations (both WDM and CDM) over a wide range of scales including very small scales up to 5 
kpc. 

It is easy to introduce the cosmological constant in the framework and equations presented here. Moreover, baryons 
and photons can be treated in this framework at the price of introducing further coupled Volterra integral equations. 

In section|ll]we derive the linearized and coUisionless Boltzmann-Vlasov (B-V) equations for DM and for neutrinos. 
We present the linearized Einstein equations for the gravitational potentials which are coupled to the B-V equations. 

In section Hill we then provide the adiabatic initial conditions for the fluctuations which turns to be constrained by 
the B-V and linearized Einstein equations. 

In section IIVI we recast the linearized DM and neutrino B-V equations as a system of linear integral equations of 
Volterra type. 

We derive in Appendix |X] the Poisson equation from the explicit solution of the linearized Einstein equations and 
the systematic corrections to it in the short wavelength regime {^dm a y ^ 1). Some useful angular integrals are 
computed in Appendix |B] 



II. THE BOLTZMANN-VLASOV EQUATION IN THE FRW UNIVERSE 

We derive in this section the coUisionless B-V equation for DM particles which decoupled being ultrarelativistic 
(UR) and become non-relativistic in the radiation dominated era. This treatment is appropriate for keV scale DM 
particles which become non-relativistic by z ~ 2 x 10^ and applies also to any DM particle candidate decoupling being 
UR during the RD era. 
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A. Particle propagation in the FRW universe including fluctuations 



We consider spatially flat FRW spacetimes with adiabatic perturbations of the metric in the conformal gauge. In 
conformal time r] the metric takes the form 



where is the Newtonian potential and (f> corresponds to the perturbation of the spatial curvature. 
The particle propagation equations follow from the Lagrangean 

1 dx" dx^ 

2 ^ ^ ' 

where x" are the contravariant particle coordinates and A is the afhne parameter on the trajectory. 
The covariant canonical momentum follows from ea. (|2.2p as 



(2.1) 



(2.2) 



dC 



Pa 



and the equations of motion take the form 



^[-dx 



dxl^ c 



(2.3) 



dPa 1 

-dX^-2P'P'^- 
The equations of motion have to be supplemented by the on-shell condition 



= -ga [3 p'^ 



(2.4) 



(2.5) 



where m is the mass of the DM particle. 

The derivative with respect to A is related to the derivative with respect to the conformal time using eq. (|2.3p for 
a = 0: 



d _ 1 d 
drj dX 



The equations of motion (|2.4I) are then 



dP] 
di] 



{P^d,4> + Pld,^ 



(2.6) 



(2.7) 



To first order in the fluctuations, it is convenient to define the momentum qj and the energy variable E{r], q) as in 
refs. UM, 



q-j ^ + pj , £;(/?, q) = ^Jvi? a?'{7]) + q^ where = , q^ = qj qj . 
The on-shell condition ea. (l2.5p becomes to first order in the fiuctuations, 



(2.8) 



We therefore have to first order in </>, 



or = pI [I - 2 ip 



1 



a?{rj) 



1-^ E{7^,q) , 



(2.9) 



E{rj,q)^-po (1-0) 



, P ^PiP 



Pt+Po^ E\rj, q)+q- + 2<j)m^ a\r]) 
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In terms of qi the equations of motion (j2.7l) to first order in the fluctuations (j) take the form 

dqi d4> f 7 , 1^ 



'^^=..#-f£;a.^ + Ca.0) (2.10) 



and therefore 



where 



dq q, dq, ^cj) / - q"^ . 

— ^ — — ^ q [ E diTp + — n, ^^(|) ] , (2.11) 

dr] q dr] atj \ h 



qi dji^ 
rii = — , rii — (5'-' rii — 1 , rii —— — . (2.12) 
q dr] 



The total derivative d(j)/ drj can be expressed in terms of partial derivatives as 

dcj) d(j) q. 



d, - 8, ' E (^-^3) 



where from eas. (|2.3p . (|2.6I) and (12. 8|) to zeroth order in the fluctuations, 

dx' _ _ 

Combining eas. (|2.1ip and (|2.13p yields 



(2.14) 



^ = - E dii; ~ ^ {S^j - Hi n-j) djcj) (2-15) 
drj orj E 

p.=q^-En.d,^. (2.16) 
drj or] 

B. The zeroth-order WDM distribution and the space-time in the RD and MD eras 

We work in the universe where radiation and dark matter are both present. The radiation and DM densities are 
given in general at zeroth order by 



Pria) = ^ , Pdmia) - gdm / E{r],q) f^"^{q) , (2.17) 



Pr_ , _ Ndm f d^q 

a* a* J (2Tr)'^ 

where pr = 0.r Pc stands for the radiation energy density today, /o™(9) is the homogeneous and isotropic zeroth order 
distribution that freezed out at decoupling, normalized as 

dq f^^^iq) = 1 • (2.18) 



gdm is the number of internal degrees of freedom of the DM particle, typically 1 < gdm ^ 4 and the normalization 
factor Ndm reproduces the DM average density today pdm — ^dm Pc as 

Ndm m gdm f 7^ /o*"(g) = ^dm Pc , hence, Ndm = ^^1^ , (2.19) 
J (27r)-* gdm m 

where ^Idm — 0.233 is the DM fraction and pc is the critical density of the Universe 

Pc = 3 Mj>i Hi = (2.518 meV)* , 1 meV = 10"^ eV , i?o = 1-5028 lO"'^^ GeV , M|,; = . (2.20) 

8 TT G 

We consider generic freezed out distribution functions fg^iq). We call Td the scale of the average momentum q at 
the zeroth order freezed-out homogeneous and isotropic distribution at decoupling. When the DM particles decouple 
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at thermal equiUbrium, Td is just the (covariant) decouphng temperature. Td is related by entropy conservation to 
the CMB temperature today and to the effective number of UR degrees of freedom at decoupling gd as 

Td = (^^^ ' Tcmfc , where T,^b = 0.2348 meV . (2.21) 

In case the decoupling happens out of thermal equilibrium, Td gives the (covariant) momentum scale of the DM 
particles at decoupling. We thus introduce the dimensionless momentum both for in and out of equilibrium decoupling, 

Q^Yd' ^^-^^^ 

which typical values are of order one. 

We now consider the dimensionless zeroth order freezed-out density /o'"(Q) and the dimensionless normalization 
constant Ndm 

/o'™(Q)=Ti/o*"(9) , rQ'dQC\Q) = l , Nd^^^^ ^^'^t^s , (2.23) 

Jo ^ d 9dm J ^ 

where we used eas. (|2.18|) and (|2.19p . For example, we have for DM fermions decoupling ultrarelativistically at thermal 
equilibrium 

fr^Q) - 3^3) ' ^'-''^ 

where C(3) — 1.2020569 Out of equilibrium freezed-out distribution functions for sterile neutrinos 

are considered in the accompanying paper (ssj . 

Eq. (|2.23p and the value of the average DM density pdm ea. (|2.20p imposes on the parameters of the DM particle: 

gdm Ndm m = 0.6988 keV ^ . 

This relation suggests that DM decoupling ultrarelativistically can have its mass in the keV scale. Moreover, an 
increasing body of evidence from the combination of theory and astronomical observations points towards DM particles 
with mass in the keV scale [l^ - [l^ : we thus take 1 keV as the reference scale for the mass of DM particles. We 
consider gd = 100 as reference value for the number gd of ultrarelativistic degrees of freedom at decoupling in thermal 
equilibrium. This corresponds to a physical decoupling temperature Tdphys = {zd + 1) Td ^ 100 GeV, Td being the 
covariant decoupling temperature. 

The normalized momenta for fermions in thermal equilibrium and for out of equilibrium sterile neutrinos are 
defined as 

/"OO /"OO 

/*" ^ / Q" C\Q) dQ , I,-: s / Q" f^{Q) dQ . (2.25) 
Jo Jo 

Explicit expressions for them are given in the accompanying paper [35] . 

From now on we use for the dimensionless one-particle energy [see (|2.8p ]. 

e{y, Q) = ^^H^ = \J [£,dmY y'^ + where a = Ueq y and S,dm = "l"*"^ ■ (2.26) 

J^d J^d 

We find from eqs.^^TBl, (|2:23l) and (12:26)) . 

= ^ 4900 — (^Y = 5520 f —) ~" (gdm Ndm)^ ■ (2.27) 

^ Td keV VlOO/ VkeV/ ^ ^ ^ 

That is, S,dm will normally be a large number ^dm ~ 5000. The parameter S.dm is the ratio between the DM particle 
mass m and the physical decoupling temperature at equilibration redshift z^q + 1 — l/a^q ~ 3200. Therefore, ^dm is 
a large number provided the DM is non-relativistic at equilibration. 

It is convenient to use the dimensionless wavenumbers 12611 



K=kr]* and a = - — = tj where 77* = . / -p- — = 143 Mpc . (2.28) 

?dm Hq m ^Ueq ildm V -"0 



2 2 Td , , ^ _ / fteq 1 
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Some useful formulas 


Ho = ^-^^ Pc , Mpi = -— , Pdm = ^dm Pc , Pr = Pc , — = — 3200 


= \ 7^ -TT - Mpc , K = kr)* = \ idm a , gdm. Ndm = 2 tt^ 

V s^i\/ -Ho / -'d 


a 

, y = — 


3200 

~ z+l 


^dm= =4900 j^^^ (^qJ = 5520(^j^^^j (gdmNam)^ , 3^1^ =2 


1 2Td AttGv*^ gdruNdruT^ 1 47r G 7?*^ ff. iV^j/) (Trf^)^ 


R^y) 








i{y, Q) = r . , hiy, y') = lQ [Kv, Q) - 


-Ky'.Q)] 




poo poo 

Itr = / Q" /o*"(Q) dQ , Q" f^{Q) , I^"" = 
Jo Jo 


= 1 




Aamiv, k) = Adm{y, a) 9<imNd^n ^^q^ ^ ^^^^^ ^ ^.(^^ g- N'^iyHTUr ^^^^ 

1 r 1 fi \ n jdm 

My,oc) = -— —Kam{y,a) + ^^K{y,a) , = + i?.(0) ~ fi.(O) = 0.727 , A(0,a) = l 


^(77,fe)=V(0,fe)^(?/,Q) , ilj{v,k)=ilj{0,k)i){y,a) , i){0,a) = l , 0(0) ~ 
a{y,a) = adm{y,a) + a^{y,a) , a{ri,k) = 'ip{0,k) a{y,a) , adm{ri,k) = ^{0,k) adm.{y,a] 


1 + I ii,(0) = 1.291 

, a^{ri,k) = ilj{Q,k) a^{y,a) 


r{y,y') = 2 {^l + y-^l + y') , s{y) = -ArgSinh (^—^ 









TABLE I: Some useful formulas. 
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Universe Event 


redshift z 


_ a + 1 3200 

^ ae, 2+1 2+1 


DM decoupling 


rp 1 

-J - 1 fi in^^ '^^ { 3d\?. 

z-d r-u ru j 


yd =i 2 X 10-12 


neutrino decoupling 


^ 6 X 10" 


0.5 X 10-« 


DM particles transition from UR to NR 


^ 1-6 X 10^ ^ 


2/tr-an. = 0.0002 

10"^ <y< 0.01 


Transition from the RD to the MD era 


==: 3200 


= 1 


The lightest neutrino becomes NR 


_ at; 

0.05 eV 


0.05 eV 

Vtrans 


Today 


20 = 


yo ~ 3200 



TABLE II: Main events in the DM, neutrinos and universe evolution. 



The free-streaming length is given by |13l . |2£ 



Hq m ^ aeq ^drn ^dm 

where the momenta if/^ are defined by eq. p.25p and therefore, 
k I 



a 



, , keV /lOO 
■^il= and =57.2kpc 



/keV^ 

50.8 kpc ) {gdmNdm)~i 

\ m 



The DM energy density is given in general by eq. (12. 17^ that we can write as 

Pdm Tloiy) 

y is defined in ea. (|2.26l) . a^q = il^/^^Af — 1/3200 is the scale factor at equilibration, 

Pdmiy) 



7^o(y) 



V ^dm 



a^y) 



and 



(2.29) 



(2.30) 



(2.31) 



(2.32) 



T^oiy) 



jdm 

i^[i + o{edmy')] 



Td 

-l-A 



o 



1 



2Sdm y V^dm y'^ 



idm y < 1 , 

, £.d7ny>^ 



(2.33) 



When ^£;„i y > 1 and the WDM particles are nonrelativistic the WDM density from eas. (|2.31l) and (12. 33^ dilutes as 
as expected. For S^dm y the WDM particles are ultra-relativistic and from eas. (|2.3ip and (|2.33p the WDM 
density dilutes as radiation as 1/a^. Eq. (j2.33l) shows that pdm{y) and Prad{y) become equal at equilibration y ~ 1 
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(up to l/^^„ corrections), as it must be. In fig. [T]we plot log^p "T^o (j/) vs. log^Q y for fermions decoupling in thermal 
equilibrium and for sterile neutrinos decoupling out of thermal equilibrium in the x model where sterile neutrinos 
are produced by the decay of a real scalar [sil, [s^]. (These particle models are analogous to those in refs. [Ill which 
consider a complex scalar field.) 



-1.5 





loglo7^o(y)' vs. log 


y for a Fermi-Dirac distribution 


- 


loglo7^o(l/) vs. logio 


y for sterile neutrinos in the x model 



















-6 -5 -4 -3 -2 -1 

FIG. 1: logj^Q TZo(y) defined in ea. (|2.32p vs. logjQ y for fermions in thermal equilibrium and for sterile neutrinos out of thermal 
equilibrium in the x model. Both freeze-out distributions give the same TZo{y) values for ^dm y ^ 5 [as in ea. (|2.33p ] while TZo{y) 
does depend on the details of the freeze-out distribution for f^m 2/ ^ 5. For 2/^1, T^o{y) takes the constant value given 
analytically in ea. (|2.33|l . 



We find from cas. (|2.3ip . (|2.32p and (|2.33p that WDM gives a small contribution of the order l/£,dm to the radiation 
density for y < 1: 

?d™ a<i /f" 

Pdmiy) = Pr{y) 'R-oiv) =^ 7 — Pr{y) ■ 



dm 



That is, the quantity gives the order of magnitude of the ratio of densities Pr{y) / Pdm{y) for 2/^1 while the 
WDM is still relativistic. 

Taking into account eq. (|2.3ip the Friedmann equation takes the form 



with the explicit solution 



< (^) =H',^r [l+7^o(^/)] , (2.34) 



When the WDM particles are UR, we find that they give small corrections of the order 1/^dm to the scale factor a{rj) 



HV) = — r\/l + T— ^ ' «W) = — • (2.36) 



Hence, eq. (|2.36p indicates a little slow down of the order l/(,dm in the expansion of the universe when the WDM 
becomes non-relativistic around 2/ — 1- When the WDM particles are NR {^dm y ^ 1) the WDM corrections are 
even smaller, of the order 1/Cdm- 
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In summary, up to small l/£,drn or 1/C^„ corrections for ^dm J/ ^ 1 or ^dm J/ ^ 1: respectively, the scale factor thus 
results from ea. (|2.35p . 



a{i]) ^ Geq y{r]) , y(ry) = ^ ^1 + ^ j , 77 = 2 77* ( ^TTy - l) 



(2.37) 



The scale factor eas. (|2.37p has the radiation dominated behavior for t] ^ rj* and the matter dominated behavior for 
r] ^ rj*. Notice that yeq = 1 and y today — 3200. 



We have for the ratio Tlo{y) 



1 + 



(2.38) 



Therefore, we can always approximate l + 7?.o(j/) by 1 + y because in the case ^dm y T^oiv) ^ 1- We will therefore 
replace 1 + TZoiy) by 1 + y in most cases. 



We obtain for h{ri) defined as 



a drj 7]* y rj* y 



l + O 



(2.39) 



where we used eq. (|2.35p and (|2.38p . 

Modes reenter the horizon when their physical wavenumber kreenter/ o. is equal to the inverse of the Hubble radius 
H — h/a, that is. 



T y 



y 1.43 IQS kpc 



(2.40) 



C. The linear and collisionless Boltzmann-Vlasov equation for DM and neutrinos 

The distribution function fdmiVT Qi ^) of the DM particles after their decoupling is described by the collisionless 
B-V equation. The distribution function is thus a constant over the particle trajectories (Liouville): 







dfdrn dfdm , dqi dfdm , dx^ df, 



dm 



drj dr] drj dqi drj dx^ 

77, Qi — Qi, 5? = being the independent variables in the distribution function. 

To linear order in the fluctuations the distribution function of the decoupled particles can be written as 

UniV, X) = Ndrn 9d,n ffil) + 9' ^) = ^dm /o"(9) 9drn 1 + *<im(?7, 9, x) , 

Terms of order higher than one in are neglected in the linear B-V equation. We have from eq. (|2.42p 

/f"(77,g,f) = Ndm 9dm ft\q) ^dm{v,q,s) • 



(2.41) 



(2.42) 



(2.43) 



Since dqi/drj, dfdm/dx^ and dfdm/dni are of order one [see eas. (|2.15p and (|2.42p ]. we can write eq. (|2.41l) to the first 
order as 



dfdm , dq dfdm , qi df 



dm 



= 



drj drj dq E dx^ 
where we used ea. (l2.14p . Inserting the linearized distribution function eq. (|2.42l) into eq. (|2.44l) yields. 



dm 



d-^dm a , dln/o" 

— h — d,-^dm + 

or] hj d hi q 



d(j) E 
df] q 



Ui di ^ 



0, 



(2.44) 



(2.45) 



13 



where we used eq. (j2.16|) . Fourier transforming 



4>ir],x) 



(2^ 



(2 7r)3 

i k-x 



(2^ 



Ak-x j'dm 



d-^k 
(2^ 



ip(ri,k) , 



(2.46) 



ft'^ir^, g, k) = Ndm 9dm ft\q) *<i™(»7, g, k) 
equation p.45p becomes [3] 



dv E ' ding 



dm 



i E 



drj 



or, equivalently 



dm 



i q- k 

d-q E{r],q) 



dln/o*" 
ding 



z E(ri,q) V 

OTy g^ 



0, 



= 



(2.47) 



(2.48) 



(2.49) 



Neutrinos are described by a distribution function /[/(yy, g*, x) obeying after decoupling a Boltzmann-Vlasov equation 
similar to ea. (|2.4ip . Neutrinos decouple in thermal equilibrium [ll, [3 redshift ~ 6 x 10^. is the comoving 
decoupling temperature of the neutrinos ~ (1 MeV/z^) ~ 0.17 10"'^ eV. 

We can linearize the B-V equation around the equilibrium zeroth order neutrino distribution as 

(77, x) = (r,) (g) + (ry, q, x) = (r,) /q^ (g) [l + (77, q, f)] , (2.50) 

where 

foiq) 



1 



3C(3) (T-)3 e?/^^+l 



(2.51) 



The normalization of the neutrino distribution eq. (|2.50l) is fixed by the neutrino energy density being a fraction Rv{ii) 
of the radiation energy density pr = Pc in the radiation dominated era 



which gives using eq. p.5ip 



^ iV.(,7) = i?.(^) , 



(2.52) 



(2.53) 



where = 7C(4)/[2C(3)] for the Fermi-Dirac distribution. The neutrino fraction Ri,(r\) changes at the temperature 
of electron-positron annhiliation (see ref. 5] and the accompanying paper [35j ) and becomes negligible in the matter 
dominated era. 

Neutrinos can be considered massless and otherwise can be neglected. Therefore (77, g, z) obeys the massless 
version of eq. p.48p 



drj 



ding 



dri 



i Ui k^ Ip 







(2.54) 



D. The linearized Einstein equations for the gravitational potentials. 



The Einstein equations for the FRW metric plus fluctuations ea. (|2.ip give for the gravitational potential at linear 
order [l|, i] 



3 h{r]) -- + fc^ 0(77, fc) + 3 h'^irj) 7/^(77, fc) = 4 tt G 0^(77) 6T°{r], k) , 
orj 



(2.55) 
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(j){ri,k) - ^p{r],k) 



= 47r G 



(2.56) 



where h{rf) is defined in ea. (l2.39p . (5Tq contains the contributions to the energy density from the photons, neutrinos 
and DM fluctuations and S(?7, k) is the anisotropic stress perturbation. 

During the RD era radiation dominates over matter and therefore the DM fluctuations are much smaller than 
the radiation fluctuations. Thus, the gravitational potential is dominated by the radiation fluctuations (photons and 
neutrinos) . The photons can be described in the hydrodynamical approximation (their anisotropic stress is negligible) . 

The tight coupling of the photons to the electron/protons in the plasma suppresses before recombination all photon 
multipoles except Qq and Oi. (The O; stem from the Legendre polynomial expansion of the photon temperature 
fluctuations 0(7], q,k) [l|). 

8o and Oi obey the hydrodynamical equations [l| 



^ + ^e,(,,fc) = ^ , 

drj arj 

This is a good approximation for the purposes of following the DM evolution f]\ . 

The energy-momentum fluctuations are the sum of the DM, photons and neutrino contributions 

xrrot M Adm{il,k) + A^{i],k) . J} ( \ ( \ n, ( m 

while only DM and neutrinos contribute to the anisotropic stress E(?7, k) 

11(77, k) = ^dmiVi k) + T.^{ri, k) . 



(2.57) 
(2.58) 



(2.59) 



Rjiri) stands for the photon fraction of the radiation and Pr(fy) for the radiation density (neutrinos plus photons). 
R-y{ri) vanishes in the MD era. 

The DM contribution to the energy-momentum tensor and to the anisotropic stress take the form 



AdmiVik) 



(2 7r)3 



E{ri,q) fr{v,q,k)^Ndrn gdn 



dh 



(2 7r)3 



E{j^,q) fS"\q)^dm{ri.q,k) , (2.60) 



'^drniV^k) 



d^q q^ 
(2^ i?(?7,g) 



1 - 3 (fc • q) 



fr(ji.q,k) 



(2.61) 



-2 Ndm gdn 



(2^)3 i?(77,g) 



(2.62) 



where P2{x) = (3x^ — l)/2 is the Legendre polynomial of order two and Adm('7, fc) stands for the DM density 
fluctuations in general (whatever ultrarelativistic, non-relativistic or intermediate regimes). 



Similarly, the neutrino contributions take the form 



A,(7y,fc) = 7V,(7y) 

I],(r/,fc) = -2N,ir)) 
The gravitational potentials (('{i]), ip{r]) thus obey 



dh 



(2 7r)3 
d^q 



(2^) 



3 qfoiq)P2{k-q) ^4v,q,k) 



3 h(v) -7^ + k^ (/)(?7, ^) + 3 h^iv) iKv, k) = -A^ G 
or] 



Adr,Av,k) + A^{r),k) 
0^(77) 



+ Aa^{Tj) p^{ri) 80(7?, fc) 



(2.63) 
(2.64) 

, (2.65) 
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C7{r], k) = (j){r], k) ~ ip{r], k) = ^^^2^) fc) + Sj.(?7, fc) where p^{7j) = R^{r]) prirj) , (2.66) 

as follows from eas. (P351) - (P?Ml) . 

In the radiation/matter domination eras the gravitational potential equation (|2.65p takes in the dimensionless 
variables y and k the form, 

y[l+TZo{y)]^ + l{Kyf ^(y, k) + [1 + 7^o(2/)] i;{y,K) = - '^''fj^ [Ad„,{y,K) + A,{y,H)] - 2 R^{y) Qo{y,K) , 



ay 3 3 a 

where k is defined in ea. (|2.28p and we used 



(2.67) 



16 n Ga\v)p,{v)= 2 R,{v)^ , 7^o(y) - . (2.68) 

V Pr[y) 



A(im(»7, fc) is connected to the customary DM density contrast S{r],k) by 0] 

+ "p.™ae,(y + l) • ^'"^^^ 
In the short wavelength limit ^ h^, eq. (|2.65p becomes the Poisson equation, as expected 

; 2 i / lA non-rclativistic . „ ^\V) + '^eg r/ fx 

k (pdmiv^k) = 4 7rGpdm o(>?, «) ■ (2.70) 

0^(77) 

In Appendix |X] we provide the explicit integral representation (jAll) to the solution of the first order differential 
equation (|2.67l) . Then, we derive the asymptotic expansion of eas. (|2.67p and (jAl[) in the k y ^ 1 (short wavelength) 
regime. We obtain in this way the Poisson equation eq. (|2.70p plus the next to leading terms in this regime in eq. (jA2p . 
Notice that the anisotropic stress a{y, k) vanishes for k j/ 3> 1. 



III. INITIAL CONDITIONS FOR THE LINEARIZED BOLTZMANN-VLASOV AND EINSTEIN 

EQUATIONS 

We investigate in this section the initial conditions for the DM linearized distribution function \I'(?7,(f, fc) solution 
of eq. (|2.49p and the gravitational potentials (t>{y, k) and ipiy, k) which obey the linearized Einstein equations (|2.65l) . 

Strictly speaking we should take the initial conditions when both neutrinos and dark matter are decoupled, namely, 
at y = 0.5 10~^ (see Table II) instead of ?/ = 0. However, setting the initial conditions at j/ = as we do here 
introduces at most an error of the order 10~^, that we can safely ignore, because both the distribution function and 
its adiabatic fluctuations (including the gravitational potentials) are regular at y = 0. 

Eq. (|2.65p yields in the 77 = limit 

V'(0,fc) ^- ^"l^/ [Ad™(0,^) + A,(0,fc)j -2i?^(0) eo(0,A) . (3.1) 

In order 4>{'q, k) and ipi^l^ k) be regular at 77 = 0, eq. (|2.66p implies that 

Sdm(0,fc) = , E,(0,fc)=0 and _^(0,^) = , — il(0,fc) = 0. (3.2) 

orj or] 

These two conditions are fullfilled provided the integrals over the directions q of 
^'dm(0, g, fc), 95'£i,„(0, q, fc)/977, 4'i/(0,<f, fc) and d'^y{Q,q,k)/dri times the Legendre polynomial P2 [k ■ qj van- 
ish in eas. (|2.62p and (|2.64p . respectively. 

In the 77—5-0 limit all fluctuation modes become superhorizon and therefore adiabatic modes must become q- 
independent except for the proportionality to the zeroth order distributions [5|. In any case, ^'dm(0, q, k) and ^1^(0, q, k) 
must be independent of the direction of q: 

^'d„(0,g,fc) =*dm(0,(Z,fc) and ^-,(0, g, fc) = *,(0, g, fc) . (3.3) 
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The linearized Boltzmann-Vlasov equation ea. (|2.49p yields to the order rf: 
i q ■ k 



M',™(0,g,fc)-^-^^(0,fc) 



din (7 



drj 



dXnq drj 



(3.4) 



and a similar expression for the neutrino distribution function. The superhorizon arguments above and eas. p.3p - p.4p 
suggest an expansion in powers of rj and i q ■ k rj for the distribution function: 



dm 



(0, q, k) = Edmiq, k) i q ■ k + Fdm{q, k) 



arj or] 
Ea. (l3.4p determines the coefficients Edm{q,k) and Fdm{q,k) as 



(0, g, fc) = £;^(g, k)iq-k + F^{q, k) . 



EdUq,k) = ^^^/^mk)-^dAO,q,k) , i^,,„(g,fc) = ^^M!!|^(0,fc). 
d In g amq orj 

Similar equations hold for Ei^{q, k) and F^{q, k). 

Eas. p.3p and p. 51) together with the integrals eas. (|2.62l) and (I2.64p guarantee that eas. p.2l) are fulfilled. 
To the first order in t] we obtain from ea. p.49p 



(3.5) 



(3.6) 



dm 



dif- 



{Q,q,k) = (iq-k)' 



'^drniO,q,k) 



—71 — VaO, k) 

amq 



dm 



+i q-k 



din g 



d In g 9?/' 



■(0,fc) 



and an analogous formula for the neutrino distribution function ^ 



(3.7) 



The knowledge of the second derivative of the distribution functions with respect to ry at 7/ = is necessary in order 
to compute the initial anisotropic stress and the difference between 0(0, k) and ^p{Q,k) from eq. (|2.66l) . 

We compute the initial DM and neutrino density fluctuations from eas. p.60p and p.63p . respectively 



A<im(0,fc) 



Ndm Qdn 

27r2 



rf'Zg'/o (9)*dm(0,<Z,fc) , A,(0,fc) 



27r2 



d<zgVo(g)*.(0,g,fc). (3., 



Inserting this result in the linearized Einstein equations p.ip at 77 = gives 



V'(o,fc) 



4 7rG?7* 



Ndm Qdn 
27r2 



dq q' foil) *dm(0,g,fc) + 



27r2 



dq q' Kiq) *.(0,g,fc) 



-2i?..,(0) 80(0, fc) 



We compute the initial value of the anisotropic stress taking the 77 — > limit in eq. p.62p with the help of eq. p.7l 



(3.9) 



^dmiV: k) _ Ndrn Qdr. 



lim 

n-^o k rj 

dm 9dm 



{2n) 



'^''^ qP2{k-q) (q-k)' 



(2 7r) 



*dm(0,(?,fc) - 



d\nf^"^ 
ding 



{0,q,k) 



^{0,k) 



frii) 



These integrals can be evaluated using eq. p.81) and 

dn{q) 



4. (^^r^2(^-^)-^^ 



with the final result 



(3.10) 



T,dm{ri,k) 2 
hm — -^-^ — = — 
ij-^o k-' r]^ 15 



Ad™(0,fc) + ^ 7Vd„grf„ V(0,fc) / dqq;' Ciq) 



Y^u{-n,k) 

hm 

r)->-0 K r] 
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A,(0,fc)+4i?,(0) Pc^{0,k) 



(3.11) 
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Inserting this result in eq. (|2.66l) gives the difference between the two gravitational potentials at the initial time 



a(0, k) = 0(0, k) - 1^(0, k) = — AdmiO, k) + A,(0, k) 

5pr L 



Ndm gd-n 
2 TT^ Pr 



dq <z'/o'™(9) 



^(0,^) (3.12) 



where we used eqs. (|237|) . ((2?52]) . ((233)) and 



1 

2 Pr 



We see from eqs. p.ll) . (|3.9p and p.l2p that all dependence on k in the initial values of 
Adr„(0, fc), Ai,(0, A;), 6o(0,fc), 5'dm(0,9, fc), ^',^(0, g, fc) and (t(0, fc) can be taken proportional to 0(0, A;). We 
can therefore factor out 7/;(0, k) from these initial values as 



*dm(O,g,fc) = 0(O,fc)cL(g) , ^'.(0,g,^) = i^(0,^)c°(g) 



(3.13) 



More generally, because the linear fluctuations evolve on an homogeneous and isotropic cosmology, the linear evolution 
equations only depend on the modulus k (as we shall see explicitly in the next section), the dependence on the k 

directions keeps factorized for all times 77. This is true for the distribution functions '^dmiv^Qj ^) ^^'^ ^lyiv^Q^ k) and 
for both gravitational potentials and (p. 

Notice that from eq.^^ 

^dm rp Pdm 

ddm ~Z 9" d — Pd7n — ^eq ~Z ; 

2 TT^ TTl 



£,dn 



and its neutrino counterpart ea. (|2.53p . 

The initial gravitational potential ilj{0,k) is a Gaussian random field with variance given by the primordial infla- 
tionary fluctuations [l|, 0, Hj] 



<V^(0,^)^(0,fc') >=^^S{k + k') 



where we can use, 
2 7r2 



fc3 



A^(fc) 



iTT^ |AoP 

9 fc3 



, A^(fc) = - A7j(fc) , A2,(fc) |Ao|^ 



(3.14) 



(3.15) 



The subscripts ^ and iz refer to the gravitational field and the scalar curvature, respectively. |Ao| stands for the 
primordial power amplitude, is the spectral index, and ko is the pivot wavenumber [27l [28|: 



|Ao| ~ 4.94 10-5 , 0.964 , fco = 2 Gpc" 

The initial value of the gravitational potential "0(0, k) can therefore be written as 



V>(o,fc) 



3^ki \ko 

where (/(fc) is a Gaussian random field with unit variance 

<g{k) g*{k') >^6{k~k') . 



9ik) 



(3.16) 



(3.17) 



A. Physical magnitudes in dimensionless variables 

From the analysis in the previous subsection we see that it is convenient to define dimensionless density fluctuations 
and dimensionless anisotropic stress fluctuations factoring out the initial gravitational potential '0(0, fc) in order to 
obtain quantities independent of the fc direction: 

Admiy,^) ^ hdm{y,a) g-'-ft V(0,fc) , A.(y,^) = A.(y,a)^^i4^MZI^(0,^), 
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(l){y,K) = ip{Q, k) (f>{y, a) , i){y, k) ^ tp{0,k) ip{y,a) and tlj{0,a)^l, 

(Jd7n{y,K) ^ tp{Q,k) adm{y,a) , (j''{y,K) = tp{0,k) a''{y,a) , a{y,K) = ip{0,k) a{y,a) 

a{y,a) = (f>{y,a) -'4>{y,a) , o-(0, a) = (^(0, a) - 1 , 60(2/, k) = ■(/'(O, fc) 60(2/, a) • 

For ultrarelativistic neutrinos in dimensionless variables we have [see eq. (I2.26p ] : 

EiTi,q)^q = Ti[ Q , e{y,Q)^Q. 

The dimensionless density fluctuations are expressed in terms of the distribution functions as 



Adm{y,K) 



17 ^ ^(0,A?) 



where we used eqs.dMD]), and dSH]). 

We find from eqs.dlSl), ([333]) and 



The customary DM and neutrino number density fluctuations are related to Ac(„i(y, k) and A,y(y, k) by 

-C'dm(2/, n) = Adm(y, k) , Nu{y, k) = A^(y, k) . 

The linearized Einstein equations (|2.67p become for the dimensionless quantities eg. p.lSp . 
d 

; + 1 ) + ^ (k y)" (P[y,a) = [l+Ko[y)\ (J[y,a)~^^-^ Adm[y,a) — — 



(l+^o(y)) ( ^ + 1 



(3.18) 
(3.19) 

(3.20) 



/ /o'"(Q) cL(Q) , A40,k)= Q' dQ f^{Q)c^,{Q) . (3.21) 

Jo Jo 



(3.22) 



Admiy,a)- ^^^l^ A^{y,a)~2R^{y) Qo{y,a) . 



(3.23) 

From eqs. p.l2p and p.lSp the dimensionless density fluctuations and anisotropic stress fluctuations take as initial 
values, 



a(0,a) = i 



A,™(0,a) + :^ A.(0,a) 



Eas. p.23p and p.24p suggest to introduce the quantities 



A(y,a) 



^drii J 3 



jdm 

7^ + i?.(0) 



A(0,a) = 1 



rdn 



_ Ad,n{0,a) R,{0) 

+ —nr- A^(0,a) . 
-'3 



dm 



The relation between the initial values ea. (l3.24p becomes, 

4 1 - 

5 CT(0,a) =4/4 +/^ and 0(0, a) = 1 + - + - . 

o o 

The linearized Einstein equations f|3.23p can be thus written in a more compact form 



(l + 7^o(2/)) ( + + 



0(2/, a) - [1 + 7^o(y)] a(2/, a) - ^ A(2/, a) - 2 i?^(y) 60(2/, a) 



Eg. (13.271) at 2/ = gives the relation 



jdra ^ 

l + -^ = --/4-2i?^(0) 60(0, a) 



(3.24) 



(3.25) 



(3.26) 



(3.27) 



(3.28) 
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where we used from eq. dOS)) that 7^o(0) = /|'"/C<im- 

The mitial number density fluctuations of photons Go(0,a), neutrinos Ni,{0,a) and DM Ddm{0,ci:) are customary 
set equal to each other P, 0, H, [HI, [H] which gives from eqs.^^, (13:25)) and ((3^281 



jdm 



J din 

2 ^ DdmiO: a) - 2 R^{0) 60(0, a) , 

<,dm 



and therefore 

N,{0, a) = 60(0, a) = Dd^iO, a) ^ 
It follows in addition from eqs. (l3.25p and (|3.30p that 



Ad™(0,a) = -2/| 



dm 



A,(0,a) = -2/3^ 



-2 1, 



-2RM 



(3.29) 



(3.30) 



(3.31) 



The approximation symbol ~ here indicates that DM contributions to the initial data of the order ^ 1 have 

been neglected. As is known, DM is negligible in the RD era and its contributions to the initial data relative to the 
radiation contribution are of the order 1 / S^dm ■ 

Using eq. p.3ip we can rewrite eq. p.26p as the relation between the two initial gravitational potentials 

2 



(/>(0,fc) 



1 



When corrections 1/^dm are neglected this becomes a known relation 0, @ 



0(0, fc) 



l + -i?.(0) 
5 



In summary this yields for the initial gravitational potential 

2 



0(0) = 0(0, a) 



1 



i?.(0) 



(3.32) 



Eas. (|3.21[) . p. 301) and p.3ip impose constraints on the functions c"^((5) and c1{Q) defining the initial distribution 
functions. We have to specify the initial functions c^dmiQ) ''^^'^ '^{Q) to completely define the initial data. There are 
two well motivated physical initial conditions. First, the thermal initial conditions (TIC) (or thermal perturbation) 



Td^Td 



1 + 



5T{k) 
Td 



in which case c|J„((3) and c^(Q) are proportional to d\x\ f^™ / d\\\Q and dln/p /dlnQ, respectively. Second, the 
Gilbert initial conditions (GIC) [l^ [2^ where c^dmiQ) ^^'^ are chosen to be constants. In order to fulfill 

ea. p.30p we must choose for both DM and for neutrinos 



1 dln/o^ 



dm 



for thermal initial conditions (TIC) , 



1 dln/o^ 



~'dm 



(Q) 



2 d\nQ 

—2 for Gilbert initial conditions (GIC) 
This completes the analysis of the initial conditions. 



2 d\nQ 
-2 for GIC 



for TIC , 



(3.33) 



IV. THE LINEAR BOLTZMANN-VLASOV EQUATION AS A SYSTEM OF VOLTERRA INTEGRAL 

EQUATIONS 

We recast in this section the linearized DM and neutrino B-V equations (|2.49|) and (I2.54p for 'i{y,q,k) and 
\l/y(?/, q, k), coupled with the linearized Einstein's equation, as a system of linear integral equations of Volterra type. 
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A. From the Boltzmann-Vlasov equations to the Volterra integral equations 



In the dimensionless variables eas. (|2.22|) - (|2.23[) the B-V equation (|2.49p takes the form 



dlnQ 



dy e{y,Q) 
It is convenient to set 

where a is related with k according to eqs. p.28p - (j2.30p and 



= 0. (4.1) 



(4.2) 



dy' 



e{y',Q) ^n + UoW) Jo 



dy' 



l + TZoiy')] y'^ + iQ/UmY 



(4.3) 



the one-particle energy e{y,Q) is defined by ea. (l2.26p . Notice that the free-streaming distance l{y,Q) depends on Q 
through the ratio Q/£,dm- From eq. ()2.33p and the discussion after it, we can set from now on TZo{y) — y ignoring 
inessential or corrections. (Except in sec. IV B of ref. [35]). 

Since q/E{ri, q) — Q/e{y, Q) is the velocity of the DM particle at time rj, its corresponding coordinate free-streaming 
length is given by 



Afs = q 



dif 



^Qi{y,Q) 



E{ry,q) ^dm 



Ifs 



Qi{y,Q) , 



(4.4) 



where we used eqs. ()2.37p and ()4.3p . Ifs is given by eq. ()2.29p and sets the scale of the coordinate free-streaming length 

Afs- 

Inserting ea. (|4.2p into ea. (|2.48p yields for ^i(?/, Q, k) the equation 



9*1 __rfln^^+.5.QK.,Q)/2 



dy 

Integrating on y we obtain: 



dlnQ 



d(j> ie{y,Q) 



dy \/l + y ' 



K • Q ^l^[y,K) 



*i(y,Q,K) = *i(0,g,K) - 



dy e 



d In Q 7o 



+ia-Q l{y',Q}/2 



oy' VT+Y 



Integrating the term d(j)/dy' by parts in y' and using eas. (j4.2p . we find for dm{y, Q, «): 



-ia-Q Hy,Q}/2 



+ 



rfln/o^'" 
dlnQ 



^-ra-Qliy,Q)/2 0(0, «) _ ^(y, a) 



+1 



K-Q dln/(f" r dy' +^s-q [i(y',Q)-i(a,Q)]/2 



g2 dinQ Jo yr+F 



£(2/',Q) 



£(2/',Q) 



^{y',a)-eiy',Q) a{y',a)] (4.5) 



Here we used eas. p.3p and (|3.13p for the initial value 5'(i,„(0, Q, k). Multiplying both sides of eq. (|4.5p by 
^iViQ) /(^"'(Q): integrating over Q and using eq. p.20p for the DM density fiuctuations yields, 



^dm{y, a) = a{y, a) + y ^dm bdmiy) (t){y, a) + n 



dy' 



vT+y 



= [N^{y,y')^[y',a)+K{y,y')a{y',a)\ , (4.6) 



where we factored out the initial gravitational potential "0(07 '?) from the density fluctuations according to eas. (|3.13p - 
p.lSp in order to obtain a quantity independent of the directions of a: 



a{y,a) 



dQe{y,Q) 



friQ)^dnAQ)+m 



dfS 
d\nQ 



Jo 



■Qi{y,Q) 



(4.7) 
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y Uni bdm{y) 



oo ^2 



Q' dQ 





Nciv, y') = / £(y, Q) ^ ji [a lQ{y, y')] 



dQ 



e{y',Q) 



e{y',Q) 



7 P (^YTh 

We used eas. (|Bl|) and (jB4|) . j„(a;) for < n < 3 are spherical Bessel functions [soj . 



(4.8) 
(4.9) 
(4.10) 



/Q(y,y')^ 2 Q[^(y,Q)-'(y',g)] = f 



da; 



(1 + a;) a;2 + (Q/ed„0^ 



and we used the relation 



Notice from eq. p.32p that 



3 a^^ 



2 7r2 



2 6 



dm 



(4.11) 



,^(0)~l + -i?,(0). 
5 

The kernels N^iy, y') and N^{y, y') only depend on the modulus of a and not on its direction since we consider linear 
fluctuations evolving on an homogeneous and isotropic cosmology. 

We derive now for adm{y, ol) an equation analogous to cq. (|4.6p . We first obtain from eas. (|2.62p . (|2.66p and p.l3p - 

(Eni), 



V'(0, a) crd,n{y, a) = ^ ^dmiri, k) ^ - 



We multiply eq. (|4.5p by 



k'^ y'^ J 47r e{y, Q) 



e{y,Q) 

integrate over Q and using eas. p.20l) and (I4.12p we find, 



P2 • Q) f^'^iQ) *<i™(y, Q, ■ (4.12) 



(4.13) 



6dm CTdm iy,a) = a''{y,a) + K 



dy' 



VI + 2/' 



= [U^{y,y')^{y',a) + U:^{y,y')a{y',a)\ 



where, 



a'^{y,a) 
Uciy.y')^ 
K{y,y"- 



dQ 

k'^ y'^ Jo £iy,Q) 



ft{Q)^UQ) + m 



dm 



d\nQ 



J2 



i^QKy.Q) 



(4.14) 



(4.15) 



Q^dQ df^ 



dm 



Bn'^y^Jo e{y,Q) dQ 



e{y\Q) 



Q' 

eiy',Q) 



{2 ji [a lQ{y, y')] - 3 j3 [a lQ{y. y')\} , (4.16) 



oo ^4 



Q" dQ dfg 



dm 



e{y', Q) {2 n [a Igiy, y')] ~ 3 h [a Iqiy, y')]} 



BK^y^Jo e{y,Q) dQ 
We used here eq.([Bl| and (jB4| . 

Eqs. (|4.6p and (|4.14l) form a system of Volterra equations 

- dy' 

Adm{y,a) = a{y,a) +y ^dmbdm{y) (p{y,a) + K / , , [Na{y,y') (l){y' ,a) + N^{y,y') a{y' ,a)] , 



(4.17) 
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fly' 

idm&dm{y,a)=a%y,a) + K [u^{y,y') 4,{y' ,a) + UZ{y,y') a{y' ^a)] . (4.18) 

Notice that a(y,a), a'^(y,Q!), (j){y,a)^ ^dm{y,ct), ^dm{y,cy) and a{y,o) only depend on the modulus of a and not 
on the directions of a. The dependence on the a. directions conies from the initial power spectrum 7/;(0, k) through the 
random field g{k) in eq. (j3.17l) and turns to factor out, which simplifies the resolution of the Volterra integral equations 
()4.18p . The factorization of the dependence on the a directions is possible because we consider linear fluctuations 
evolving on an homogeneous and isotropic cosmology where all the evolution kernels Na(jj,y'), N^{y,y'), Ua{y,y') 
and U^{y,y') are independent of the a directions. 

The B-V distribution function as well as the coefficients in the B-V equation depend on y, a and Q. We integrate the 
distribution function over Q multiplied by appropriated weights. The distribution function times e{y, Q) produces the 
density eq. (|3.20p and the distribution function times the expression (|4.13p produces the anisotropic stress fluctuations 
eq. (|4.12p . The density and the anisotropic stress fluctuations defined with such specific weights obey a closed system 
of Volterra integral equations. Namely, no extra information on the Q dependence of the distribution functions is 
needed, which is a truly remarkable fact. 

We derive below the Volterra integral equations for neutrinos eas. (|4.26p similar to eas. (|4.18p for DM. 



B. The pair of Volterra integral equations for DM and neutrinos 

The Volterra integral equations for neutrinos are obtained from eq. ()2.54p following the same steps eqs. ()4.2p - (j4.17p 
which lead to the DM Volterra integral equations (|4.18p . These Volterra integral equations for ultrarelativistic 
neutrinos are simpler than the corresponding DM equations and follow from eas. (|4.7p - (|4.17p making the following 
substitutions: 

£(y,g)^g , alQ{y,y')^ nr{y,y') , a Z(y, Q) ^ ^ r(y, 0) , gdm ^ {y) , (4.19) 

fo""{Q)^ foiQ) , Adrn{y,a)^A,{y,a) , , admiy,a) ^ a-'iy^a) = r{y,c^) - riy,a) , 

where 

r(y,2/') = 2 (v/TT^-yTT7) , r(y, 0) = 2 (VT+i; - l) (4.20) 

and we used eq. (|4.3l) . [See also eq. (|3.19p ]. 

Upon these changes the kernels A^Q(y, y'), N^{y,y'), Ua{y,y'), Ugiv- v'^ in eqs. (l4.9p - (|4.17p and the inhomogeneous 
terms a{y,a) and a"{y,a) in ea. (|4.7p simplify considerably. For ultrarelativistic neutrinos (ur) (0 < y < 34 mi,/0.05 
eV) using eas. p.l9l) - (|4.20p these kernels become: 

iVo(2/, y') " (y, y') ^ -8 n r(y, y')] , 

U^{y, y') "°=^~ C/r (y, y') ^ {2 Ji [n r{y, y')] - 3 js r{y, y')]} , (4.21) 

5 K 

K{y,y')^'''n''-~lNZ{y,y') , f/:(y,y')""^=^'''"°^-^C/r(2/,y'), (4-22) 

a{y, a) " ""S^'™"^ a^^iy, a) -2 [l + 2 ^(0)] jo r{y, 0)] , (4.23) 

a-(y,a)"'^'==^""°^a"'-'^(y,a)^-6/3-[l + 2^(0)] [^[^^j")] , (4.24) 

K y 



where we used eqs.(I22Sl), (E2II), (0301), and (|il^ - (|iT7p . 

2 massless and using Tab 

y(dmbd„,{y)^4I^ . (4.25) 



In addition, when relevant the neutrinos are massless and using Table U the coefficient of (j){y,a) in eqs. (|4.18p for 
neutrinos becomes: 
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Therefore, making the changes eas. (|4.19|) - (l4.25|) in eas. (|4.18|) yields the following Volterra integral equations for 
ultrarelativistic neutrinos 



Jo yi + y 



<l>{y\oi) - i (j{y',a) 



(4.26) 



-'3 



R. 



^a''{y,a)=a-^^y^a) + K -^^UZ{y,y') 

Ay) Jo VI + y 



dy' 



(4.27) 



Notice that the DM and neutrino Volterra integral equations eqs. (j4.18p and (|4.26l) are coupled to each other and to 
the linearized Einstein equations eq. (|3.23p as well as to the hydrodynamic photon equations (|2.57l) - (|2.58l) . 

It is possible to simplify the set of four Volterra integral equations (|4.18l) and (|4.26p into two Volterra equations. 
Taking linear combinations of eqs. (|4.18p and (|4.26p we find for A{y, a) [defined in eq. p.25p ] and (j{y, a), 

A(y,a) = C(y,a) +i?5(y) ^iy,a) + f dy' [Go.{y,y') ^{y' ,a) + Gl{y,y') a(y',a)] , 



a(2/,a)=C"(y,a)+ / dy' [r^{y,y') a{y' ,a) + I^{y,y') ^{y' ,0)] 
Jo 

with the initial conditions eas. p.25p and p.3ip 



A(0,a) = l , a(0,a) = ^ ~ ^ i?,(0) 



We have in ea. (|T^ 

C{y,a) = 

Ga(2/,y') = 
Gl{y,y')-- 
Ia{y,y') = 



1 

"27^ 
1 

2ic 



a(y,a) Ry{y) „ , 



[yhdm{.y) + ^R,.{v)\ , 
1 



, C%y,a) 



_ a'^iy,a) Ry{y) 



2/4 VT+F 



NM) + ^N:'{y,y') 



J^Kiy,y')-^N:^iy,y') 



R.{y) 



— Ua.{y,y') + ^K^{y,y') 



^K{y,y')~^K^{y,y') 



^d 



^3 



ur a I \ 

a (y,a) 



(4.28) 



(4.29) 

(4.30) 
(4.31) 
(4.32) 
(4.33) 



In eqs. (|4.28l) we can use ~ Riy{0)- 

Notice that the G and / kernels in eas. (|4.30l) - (|4.33l) result expressed as the sum of the DM contribution from 
the and U kernels plus the (ultrarelativistic) neutrino contribution N^'^{y,y') and U^'^{y,y'), respectively. The 
inhomogeneous terms G{y,a) and C"^(y,a) and the coefficient B^{y) in eqs. (|4.29p - (|4.30p turn also to be expressed as 
the sum of the DM plus the neutrino contributions. 

In the MD era the neutrinos are negligible and its fraction R^{y) becomes <C 1 and can be neglected. Once neutrinos 
are negligible in the MD era, the DM contribution to f7(y, a) from eqs. (|4.28p - (|4.33p is of the order ^ 1 and the 

anisotropic stress becomes negligible. This reduces the coupled Volterra integral equations (|4.28p in the MD era to a 
single Volterra integral equation for A(?/, a) as we explicitly show in the accompanying paper [35]. 

All functions in the inhomogeneous terms, coefficient and kernels in the the Volterra equations (|4.28p are explicitly 
known from eqs. (|4.29p - (|4.33p Therefore, eqs. (14.281) plus the linearized Einstein equations p.27p and the hydrodynamic 
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photon equations (j2.57p - (|2.58p provide a close system of equations determining A{y,a), 4i{y,a) and a{y,a). Once 
A(y,a), (j){y,a) and a{y,a) are known we can insert them in the r. h. s. of eas. (|4.18p and (|4.26p to obtain 
Adm(j/,a), adjn{y,a), A''{y,a) and a" {y , a) , respectively. 

We now set y = in the system of Voltcrra equations (|4.28p to check their consistency. Taking into account eqs. 
(Ii7)) - (HJU)) . (|il^ - (|iTfl) and we obtain 

C(O,a) = l + 2 0(O) , B^{0) = -2 , ^^(0, a) = [l + 2 0(0)] , 

hm^ r dy' G^[y,y')^{y\a)^Q , lim^ f dy' Gl{y,y') a{y\a) = G , (4.34) 



hm r dy' Uy,y')^{y',a)^\l^m , 1™ C dy' rjy,y') &{y\a) = . (4.35) 
a^o Jo 5 y->o _/o 25 

Eqs. (|4?28)) are identically satisfied at y = due to eas. ((3?26t and (|434)) - (|435)) . 

The system of Volterra equations (|4.28p is collisionless and it is therefore valid after both DM and neutrinos 
decoupled for j/ > ~ 0.5 10"^ (see Table II). Since we are interested in adiabatic fluctuations which are regular 
solutions of eas. (|4.28p at y = we can start the evolution at y = instead of y = ~ 0.5 10^® with a negligible 
error. 

For the DM particles, the range 0.5 10^^ < y < 0.01 corresponds to the transition from ultrarelativistic to non- 
relativistic kinematics (see Table II). 

The density contrast (5(y, a) can be expressed in terms of the normalized DM fluctuations Adm(y, a) from eqs. (|2.69p 
and (PIB as 

%,a) = -^ ^it^ 5(0,d) = -^^(0,c5), (4.36) 

where we used eq. (|3.30p and ■0(0, a) is given by the primordial fluctuations ea. p.lTp and ^dm is given explicitly by 
eg. ([2:271) . 

The integral equation (|4.28p supplemented by the fluid equations (|2.57p - (|2.58p for the photons and the linearized 
Einstein equations ()3.23p provide a closed system of equations to determine the DM, photon and neutrino density 
fluctuations. This system of Voltcrra-typc integral equations is valid for relativistic as well as non-relativistic particles 
propagating in the radiation and matter dominated eras. This is the generalization of Gilbert's equation which is 
only valid for non-relativistic particles in a matter dominated universe [19[. 

We solve in an accompanying paper [35j the cosmological evolution of warm dark matter (WDM) density fluctuations 
presented here in the absence of neutrinos. In that case the anisotropic stress vanishes and the Volterra equations 
(|4.28p reduce to a single integral equation. 
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Appendix A: The DM gravitational potential for large wavenumbers. 

In sections lU and HVl we found integrals of the type 

h{y) = f ^ {\^e--A' f{x) , X = nVi (Al) 
7o l + x \l + x J 

The function 4>\{y) = I\{y)/y solves the first order differential equation 

d 



y{l + y) ^ + i + v + \y'^ 

dy 



(t)\{y) = f{y) 

which has the form of the linearized Einstein equations (I3.23P and p.27p . 
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We derive here the asymptotic expansion of I\{y) in the hmit where A ^ 1. 
It is convenient to change the integration variable x in eq. (|Al|) to s defined as 



s{x) = loj 



The integral in eq. ()Aip becomes 



l + x 
1 + 2/ 



+ y-x , s(2/) = , s(0) = 2/ - log(l + ?/) 



ix(y) 



y-log(l+y) 7 

e-'\f[x{s)]^ 
x{s) 



In the A ^ 1 regime this integral is dominated by the end-point of integration s = 0. Expanding f[x{s)]/x{s) around 
s — and integrating term by term yields 



A>i f{y) _ 1 + 2/ 

Xy {Xyr 



df f{y) 
dy y 



O 



1 



(A yy 



(A2) 



Appendix B: Angular Integrals 

We proceeded in sec. IIVI to compute integrals over the directions of q with the help of the partial wave expansion 



1=0 

Integrating this expansion over the angles yields [SC 

dn{q) 



OO 



We use in sec. |IV] eq. dBl]) for < / < 3, the relations [2S 



e'^'^^P, (Bl) 



/ N sin a; . , ^ ^ , ^ dji , „ /t,„x 

jox = , jz+i X = - X - -f , />0 B2 

X x ax 

and the formulas for Legendre polynomials (30l | 

Poix) = 1 , Pi{x) = X , Pi+i{x) = X Pi{x) - Pi-i{x) , 1>1. 

It follows from these relations, in particular, that 

xP2{x)^lp3{x) + lPiix). (B3) 
5 5 



We get combining eas. (|Bip and (|B3I) , 



dn{Q) 



+^''-<i^'^y''Q^-'<~y''i^y'k-QP2{k-Q)=-^ {2n[alQiy,y')]~3j3[alQ{y,y')]} . (B4) 
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